# simplified version of StrategyUnitTheme::su_theme_pal()
# see: https://github.com/The-Strategy-Unit/StrategyUnitTheme
su_theme_pal <- function(palette = "main", ...) {
pal <- c(orange = "#f9bf07",
red = "#ec6555",
blue = "#5881c1")
grDevices::colorRampPalette(pal, ...)
}
If you haven’t already got these packages installed you will need the following to complete this workshop:
install.packages(c("tidyverse", "simmer", "simmer.plot", "gridExtra", "TruncatedNormal"))
suppressPackageStartupMessages({
library(tidyverse)
library(simmer)
library(simmer.plot)
library(TruncatedNormal)
library(writexl)
library(gridExtra)
})
For the patients interarrival rate, the time taken for asking the screening questions and the time for the receptionists to make the decision on what appointment to book.
You can use helper functions
# generic truncated normal distributions looks like:
rtnorm(n = 1, mu = 1, sd = 0.25, lb = 0, ub = Inf)
## [1] 0.6213012
dist_helper <- function(rfn, ...) {
function() rfn(1, ...)
}
dist_patient_interarrival <- dist_helper(rexp, 50 / 60)
dist_screening <- dist_helper(rtnorm, 4, 1, 0, Inf)
dist_decision <- dist_helper(rtnorm, 1, 0.25, 0, Inf)
dist_booking <- dist_helper(rtnorm, 1, 0.25, 0, Inf)
Alternatively:
dist_patient_interarrival <- function() rexp(1, 50 / 60)
dist_screening <- function() rtnorm(n = 1, mu = 4, sd = 1, lb = 0, ub = Inf)
dist_decision <- function() rtnorm(n = 1, mu = 1, sd = 0.25, lb = 0, ub = Inf)
dist_booking <- function() rtnorm(n = 1, mu = 1, sd = 0.25, lb = 0, ub = Inf)
dist_contact_mode <- function() sample(1:2, 1, FALSE, c(0.15, 0.85))
dist_booking_type <- function() sample(1:3, 1, FALSE, c(0.3, 0.6, 0.1))
Think about the flow diagram of the process to create the trajectories for the model, let’s break it down into manageable sections.
Create a trajectory for this branch called branch_ecr. Within the trajectory use set_attribute via_ecr to label the patient with a value of 1, so we will know what route they have come to analyse later. Then seize a receptionist ready for the next stage of the process.
branch_ecr <- trajectory("branch ecr") %>%
set_attribute("via_ecr", 1) %>%
seize("receptionist")
branch_ecr
## trajectory: branch ecr, 2 activities
## { Activity: SetAttribute | keys: [via_ecr], values: [1], global: 0, mod: N, init: 0 }
## { Activity: Seize | resource: receptionist, amount: 1 }
create a trajectory for this branch called branch_phone. Within the trajectory use set_attribute via_ecr to label the patient with a value of 0 (as opposed to 1 which are ECR route patients).
The patients renege after 5 minutes so log this and set the attribute "reneged" to 1. Then seize a receptionist ready for the next stage. Once the patient has a receptionist abort the renege and add timeout of distribution screening
renege_hungup <- trajectory("hung up") %>%
set_attribute("reneged", 1) %>%
log_("lost my patience. Hanging up...")
branch_phone <- trajectory("branch phone") %>%
set_attribute("via_ecr", 0) %>%
renege_in(t = 5, out = renege_hungup) %>%
seize("receptionist") %>%
renege_abort() %>% # now have a receptionist the patient won't renege
timeout(dist_screening)
Receptionist is deciding if a patient needs phone appointment, needs phone appointment or does not need any appointment.
branch_booking_none <- trajectory("no appointment needed") %>%
release("receptionist") %>%
set_attribute("booked_appt", 0)
branch_booking_none
## trajectory: no appointment needed, 2 activities
## { Activity: Release | resource: receptionist, amount: 1 }
## { Activity: SetAttribute | keys: [booked_appt], values: [0], global: 0, mod: N, init: 0 }
branch_booking_phone <- trajectory("booking a phone") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 1)
branch_booking_phone
## trajectory: booking a phone, 3 activities
## { Activity: Timeout | delay: function() }
## { Activity: Release | resource: receptionist, amount: 1 }
## { Activity: SetAttribute | keys: [booked_appt], values: [1], global: 0, mod: N, init: 0 }
branch_booking_f2f <- trajectory("booking a f2f") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 2)
branch_booking_f2f
## trajectory: booking a f2f, 3 activities
## { Activity: Timeout | delay: function() }
## { Activity: Release | resource: receptionist, amount: 1 }
## { Activity: SetAttribute | keys: [booked_appt], values: [2], global: 0, mod: N, init: 0 }
patient_flow <- trajectory("patient flow") %>%
branch(dist_contact_mode, TRUE,
branch_ecr,
branch_phone) %>%
timeout(dist_decision) %>%
branch(dist_booking_type, TRUE,
branch_booking_none,
branch_booking_phone,
branch_booking_f2f)
plot(patient_flow, fill = su_theme_pal("main"))
# you can change the verbose argument to TRUE to get more information
plot(patient_flow, fill = su_theme_pal("main"), verbose = TRUE)
You can use to() function in the generator, so that it completes the arrivals. Use to() with caution as the stop time is that of the generator and not the simulation. The simulation will run until all arrivals are finished, so we are assuming that staff stay overtime
env_part1 <- simmer("part 1") %>%
add_resource("receptionist", 2) %>%
add_generator("patient", patient_flow, to(8 * 60, dist_patient_interarrival), mon = 2)
env_part1 %>%
get_mon_arrivals() %>%
arrange(start_time) %>%
head(10)
| name | start_time | end_time | activity_time | finished | replication |
|---|---|---|---|---|---|
| patient0 | 0.9062182 | 8.666809 | 7.760591 | TRUE | 1 |
| patient1 | 1.0810663 | 2.880660 | 1.799594 | TRUE | 1 |
| patient2 | 1.6043486 | 7.786732 | 4.906072 | TRUE | 1 |
| patient3 | 2.0532035 | 7.053204 | 0.000000 | TRUE | 1 |
| patient4 | 2.3693078 | 7.369308 | 0.000000 | TRUE | 1 |
| patient5 | 2.7333249 | 7.733325 | 0.000000 | TRUE | 1 |
| patient6 | 3.6477607 | 13.937513 | 6.150780 | TRUE | 1 |
| patient7 | 3.9308550 | 13.606157 | 4.939348 | TRUE | 1 |
| patient8 | 4.9175231 | 9.917523 | 0.000000 | TRUE | 1 |
| patient9 | 7.7549414 | 12.754941 | 0.000000 | TRUE | 1 |
env_part1 %>%
get_mon_arrivals() %>%
arrange(start_time) %>%
tail(10)
| name | start_time | end_time | activity_time | finished | replication | |
|---|---|---|---|---|---|---|
| 389 | patient388 | 472.4713 | 477.9166 | 5.445350 | TRUE | 1 |
| 390 | patient389 | 472.8399 | 483.2193 | 6.273473 | TRUE | 1 |
| 391 | patient390 | 473.1776 | 481.5829 | 3.666322 | TRUE | 1 |
| 392 | patient391 | 474.4332 | 479.4332 | 0.000000 | TRUE | 1 |
| 393 | patient392 | 475.5632 | 480.5632 | 0.000000 | TRUE | 1 |
| 394 | patient393 | 476.6170 | 487.7046 | 6.121711 | TRUE | 1 |
| 395 | patient394 | 476.7175 | 481.7175 | 0.000000 | TRUE | 1 |
| 396 | patient395 | 477.5498 | 482.5498 | 0.000000 | TRUE | 1 |
| 397 | patient396 | 478.8235 | 490.1992 | 6.979876 | TRUE | 1 |
| 398 | patient397 | 479.2955 | 484.2955 | 0.000000 | TRUE | 1 |
If you are new to R and struggling with data manipulation, you can export the results as xls or csv. Or email us and we can help you!
write_xlsx(get_mon_arrivals(env_part1), "arrivals_part1.xlsx")
write_xlsx(get_mon_resources(env_part1), "resources_part1.xlsx")
write_xlsx(get_mon_attributes(env_part1), "attributes_part1.xlsx")
env_part1 %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
summarise(countreneged = n())
| countreneged |
|---|
| 179 |
n = 148 with 2 receptionists
env_part1 %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
summarise(countF2F = n())
| countF2F |
|---|
| 30 |
n = 27 with 2 receptionists
Get attributes of finished patients from get_mon_attributes(), join with get_mon_arrivals(). Because of varying number of decimal places we get some rounding errors resulting in near zero but negative results. That’s why when we are doing data wrangling, wee use pmax (), so it takes 0 or the wait time which ever is greater. We can then calculate summary statistics.
part1_attributes <- env_part1 %>%
get_mon_attributes() %>%
dplyr::select(-time) %>%
pivot_wider(names_from = "key", values_from = "value", values_fill = 0)
part1_callwait_results <- env_part1 %>%
get_mon_arrivals(ongoing = TRUE) %>% # not needed as there aren't any due to the generator and run()
inner_join(part1_attributes, by = c("name", "replication")) %>%
filter(via_ecr == 0) %>%
mutate(wait_time_call = case_when(
# the caller reneged
reneged == 1 ~ 5,
# in queue still - these could also still be in progress
# (would need to know time resource seized to calculate wait time),
# except we used to() in generator and run() so we have no unfinished
# !finished ~ (10 * 60) - start_time,
# everything else, (assume has been answered)
TRUE ~ pmax(end_time - start_time - activity_time, 0)
))
part1_callwait_summary <- part1_callwait_results %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(!finished),
countfinished = sum(finished))
part1_callwait_summary
| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|
| 5 | 4.223388 | 1.733617 | 0 | 5 | 3.906749 | 5 | 0 | 329 |
with 2 receptionists
# median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait
#1 4.839153 3.811972 2.657451 0 5 3.010095 5
# countunfinished countfinished
#1 0 334
Glimpse the end tail and also check whether ECR backlog is cleared at the end of the days. Overall backlog - let’s look at 1 replication
attributes <- env_part1 %>%
get_mon_attributes() %>%
filter(replication == 1)
test <- env_part1 %>%
get_mon_arrivals() %>%
filter(finished == FALSE, replication == 1) %>%
inner_join(attributes, by = "name") %>%
arrange(start_time)
test
| name | start_time | end_time | activity_time | finished | replication.x | time | key | value | replication.y |
|---|
env_part1 %>%
get_mon_attributes() %>%
filter(replication == 1, key == "via_ecr", value == 1) %>%
count()
| n |
|---|
| 69 |
You can see that everyone finished and we had total of 66 ECR.
Create a new environment
env_part2 <- simmer("part2")
env_part2
## simmer environment: part2 | now: 0 | next:
## { Monitor: in memory }
To include time dependent interarrivals create a different distribution, one for each time period (3 in total). You may also find using these helper functions to switch between minutes, hours and days useful:
t_minute <- function(t) t
t_hour <- function(t) t_minute(t) * 60
t_day <- function(t) t_hour(t) * 24
dist_patient_interarrival1 <- function() rexp(1, 120 / t_hour(2))
dist_patient_interarrival2 <- function() rexp(1, 240 / t_hour(6))
dist_patient_interarrival3 <- function() rexp(1, 40 / t_hour(2))
We are also going prioritise the patients contacting the surgery via the phone.
In the ECR branch, no changes are needed. The default priority for all arrivals is set up to 0.
In the phone branch, we want to set a higher priority using set_prioritization(), that receptionists prioritise answering calls over reviewing electronic requests but without interrupting current tasks (prioirity = 1, preemptive = 1, restart = FALSE). If you didn’t use now(env_part1) to set attributes you can use join(branch_phone) instead to link after the prioritisation back to your original phone branch.
branch_phone2 <- trajectory("branch phone prioritised") %>%
set_prioritization(c(1, 1, FALSE)) %>%
join(branch_phone)
plot(branch_phone2, fill = su_theme_pal("main"))
from the help page for add_generator():
priority: the priority of each arrival (a higher integer equals higher priority; defaults to the minimum priority, which is 0.
preemptible: if a seize occurs in a preemptive resource, this parameter establishes the minimum incoming priority that can preempt these arrivals (an arrival with a priority greater than preemptible gains the resource). In any case, preemptible must be equal or greater than priority, and thus only higher priority arrivals can trigger preemption.
restart: whether the activity must be restarted after being preempted.
same as Part 1 but with the new prioritised phone branch
patient_flow2 <- trajectory("patient flow2") %>%
branch(dist_contact_mode, TRUE,
branch_ecr,
branch_phone2) %>% #<-- changed the phone branch
timeout(dist_decision) %>%
branch(dist_booking_type, TRUE,
branch_booking_none,
branch_booking_phone,
branch_booking_f2f)
plot(patient_flow2, fill = su_theme_pal())
Use the schedule(timetable, values, period) function and set it to repeat the schedule every 24 hour period. Work out from the numbers on each shift what times capacity changes and how many are working at that point:
e.g. 0-8am = 0, 8-9am = 1, 9-9.30 am etc. becomes:
sch_receptionists <- schedule(
timetable = t_hour(c(0.0, 8.0, 9.0, 9.5, 12.0, 15.5, 18.0)),
values = c(0, 1, 2, 4, 3, 1, 0),
period = t_hour(24)
)
Add the receptionists as a resource using the schedule you created for the capacity. Either add three different generators to the model environment setting the times they operate with from_to or use functional programming to control this. You can look back at the training slides and examples to help you.
env_part2 %>%
add_resource("receptionist", capacity = sch_receptionists) %>%
add_generator("patient[8am-10am]", patient_flow2,
from_to(t_hour(8), t_hour(10), dist_patient_interarrival1, every = t_day(1)),
mon = 2) %>%
add_generator("patient[10am-4pm]", patient_flow2,
from_to(t_hour(10), t_hour(16), dist_patient_interarrival2, every = t_day(1)),
mon = 2) %>%
add_generator("patient[4pm-6pm]", patient_flow2,
from_to(t_hour(16), t_hour(18), dist_patient_interarrival3, every = t_day(1)),
mon = 2)
## simmer environment: part2 | now: 0 | next: 0
## { Monitor: in memory }
## { Resource: receptionist | monitored: TRUE | server status: 0(0) | queue status: 0(Inf) }
## { Source: patient[8am-10am] | monitored: 2 | n_generated: 0 }
## { Source: patient[10am-4pm] | monitored: 2 | n_generated: 0 }
## { Source: patient[4pm-6pm] | monitored: 2 | n_generated: 0 }
Set the random seed to 1 to give us the same results each time. Run your model for 5 days (52460) or use one of the time helper functions)
Repeat your analysis of call waits, reneging and F2F appointments as you did for part 1 using the environment for part 2. Plot utilisation of the resources using the plot function with get_mon_arrivals() and the metric “waiting_time”.
env_part2 %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
summarise(countreneged = n())
| countreneged |
|---|
| 582 |
n = 582 with 2 receptionists
env_part2 %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
summarise(countF2F = n())
| countF2F |
|---|
| 160 |
n = 160 with 2 receptionists
part2_attributes <- env_part2 %>%
get_mon_attributes() %>%
dplyr::select(-time) %>%
pivot_wider(names_from = "key", values_from = "value", values_fill = 0)
part2_callwait_results <- env_part2 %>%
get_mon_arrivals(ongoing = TRUE) %>%
full_join(part2_attributes, by = c("name", "replication")) %>% # or can use inner_join
filter(via_ecr == 0) %>%
mutate(wait_time_call = case_when(
reneged == 1 ~ 5,
TRUE ~ pmax(end_time - start_time - activity_time, 0)
))
part2_callwait_summary <- part2_callwait_results %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(!finished),
countfinished = sum(finished))
part2_callwait_summary
| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|
| 3.825883 | 3.183036 | 3.479728 | 0 | 5 | 1.464796 | 5 | 0 | 1726 |
# median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait
#1 3.825883 3.183036 3.479728 0 5 1.464796 5
# countunfinished countfinished
#1 0 1726
Some ideas:
Below are some ideas on additional analysis you might want to do
Plot waiting time over time
plot(get_mon_arrivals(env_part2), metric = "waiting_time")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
boxplot to show variation in metrics e.g. call_wait
g1 <- ggplot() +
geom_boxplot(data = part2_callwait_results,
aes(x = "rep=1", y = wait_time_call),
outlier.shape = 19,
width = 0.3,
fill = "white",
colour = "orange") +
geom_point(data = part2_callwait_summary,
aes(x = "rep=1", y = mean_wait),
shape = 19,
size = 2,
colour = "red") +
labs(x = "Part2",
y = "call wait (mins)",
title = "Boxplot of call wait from single replication") +
theme(axis.text = element_text(size = 11))
g1
Glimpse the end tail and also check whether ECR backlog is cleared at the end of the days. Overall backlog:
attributes <- env_part2 %>%
get_mon_attributes() %>%
filter(replication == 1)
test <- env_part2 %>%
get_mon_arrivals() %>%
filter(finished == FALSE, replication == 1) %>%
inner_join(attributes, by = "name") %>%
arrange(start_time)
test %>%
filter(key == "via_ecr", value != 1)
| name | start_time | end_time | activity_time | finished | replication.x | time | key | value | replication.y |
|---|
env_part2 %>%
get_mon_attributes() %>%
filter(replication == 1, key == "via_ecr", value == 1) %>%
count()
| n |
|---|
| 357 |
test %>%
head(10) %>%
as_tibble()
| name | start_time | end_time | activity_time | finished | replication.x | time | key | value | replication.y |
|---|---|---|---|---|---|---|---|---|---|
| patient[4pm-6pm]191 | 6774.405 | NA | NA | FALSE | 1 | 6774.405 | via_ecr | 1 | 1 |
| patient[4pm-6pm]192 | 6775.585 | NA | NA | FALSE | 1 | 6775.585 | via_ecr | 1 | 1 |
| patient[4pm-6pm]201 | 6790.286 | NA | NA | FALSE | 1 | 6790.286 | via_ecr | 1 | 1 |
All unfinished are ECR, because calls were prioritised ahead of them. 357 total ECR.
# A tibble: 3 x 10 - 3 are unfinished in a week out of 357
# name start_time end_time activity_time finished replication.x time key value replication.y
#<chr> <dbl> <dbl> <dbl> <lgl> <int> <dbl> <chr> <dbl> <int>
# patient[4pm-6pm]191 6774. NA NA FALSE 1 6774. via_ecr 1 1
# patient[4pm-6pm]192 6776. NA NA FALSE 1 6776. via_ecr 1 1
# patient[4pm-6pm]201 6790. NA NA FALSE 1 6790. via_ecr 1 1
Look at resources (reception staff in this example)
resources2 <- get_mon_resources(env_part2)
plot(resources2, metric = "usage", "receptionist", items = "server", steps = TRUE)
plot(resources2, metric = "utilization", "receptionist")
# show the top 10 results
head(resources2, n = 10)
| resource | time | server | queue | capacity | queue_size | system | limit | replication |
|---|---|---|---|---|---|---|---|---|
| receptionist | 480.0000 | 0 | 0 | 1 | Inf | 0 | Inf | 1 |
| receptionist | 480.0000 | 1 | 0 | 1 | Inf | 1 | Inf | 1 |
| receptionist | 481.1816 | 1 | 1 | 1 | Inf | 2 | Inf | 1 |
| receptionist | 481.6397 | 1 | 2 | 1 | Inf | 3 | Inf | 1 |
| receptionist | 484.5346 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| receptionist | 484.9699 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| receptionist | 485.1100 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| receptionist | 486.1832 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| receptionist | 486.6397 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| receptionist | 487.5739 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
# list of hours
HrList <- data.frame(
resource = "receptionist",
time = seq(from = 0, to = t_day(5), by = 60),
server = NA,
queue = NA,
capacity = NA,
queue_size = Inf,
system = NA,
limit = Inf,
replication = 1
)
# which aren't in the data
AddSnaps <- filter(HrList, !(time %in% resources2$time))
resources2_hourly <- bind_rows(resources2, AddSnaps) %>%
arrange(time) %>%
mutate(server = ifelse(is.na(server), lag(server), server),
queue = ifelse(is.na(queue), lag(queue), queue),
capacity = ifelse(is.na(capacity), lag(capacity), capacity),
system = ifelse(is.na(system), lag(system), system))
resourceutilisation_byhour <- resources2_hourly %>%
arrange(time) %>%
mutate(slot_hour = floor(time / 60)) %>%
mutate(slot_day = ceiling(slot_hour / 24)) %>%
mutate(slot = slot_hour - ((slot_day - 1) * 24)) %>%
group_by(slot) %>%
mutate(dt = time - lag(time),
capacity = ifelse(capacity < server, server, capacity),
in_use = dt * lag(server / capacity)) %>%
summarise(utilised_percent = sum(in_use, na.rm = TRUE) / sum(dt, na.rm = TRUE) * 100)
resourceutilisation_byhour
| slot | utilised_percent |
|---|---|
| 1 | 0.0000000 |
| 2 | 0.0000000 |
| 3 | 0.0000000 |
| 4 | 0.0000000 |
| 5 | 0.0000000 |
| 6 | 0.0000000 |
| 7 | 0.0000000 |
| 8 | 100.0000000 |
| 9 | 100.0000000 |
| 10 | 99.9259080 |
| 11 | 63.7850307 |
| 12 | 99.6780740 |
| 13 | 91.7546945 |
| 14 | 75.7575154 |
| 15 | 99.8493045 |
| 16 | 99.8748793 |
| 17 | 100.0000000 |
| 18 | 0.1985734 |
| 19 | 0.0000000 |
| 20 | 0.0000000 |
| 21 | 0.0000000 |
| 22 | 0.0000000 |
| 23 | 0.0000000 |
| 24 | 0.0000000 |
Option B. If you had used now(env_part1) to set attributes then you will need to copy your code and change the environment for all parts of the trajectory in which it was used. You then need to create the overall trajectory again using the new version of the phone branch (and other parts if you have used now()).
Firstly, create a new environment.
env_part2B <- simmer("part 2B")
env_part2B
## simmer environment: part 2B | now: 0 | next:
## { Monitor: in memory }
Change references to part env_part2B in all trajectories
renege_hungup2B <- trajectory("hung up") %>%
set_attribute("reneged", 1) %>%
set_attribute("reneged_at", function() now(env_part2B)) %>% # changed to env_part2B
log_("lost my patience. Hanging up...")
branch_phone2B <- trajectory("branch phone") %>%
set_prioritization(c(1, 1, FALSE)) %>% # set to increase priority from default of 0 to 1, no preempt and no restart
set_attribute("via_ecr", 0) %>%
set_attribute("arrival", function() now(env_part2B)) %>% # changed to env_part2B
renege_in(t = 5, out = renege_hungup) %>%
seize("receptionist") %>%
set_attribute("call_taken", function() now(env_part2B)) %>% # changed to env_part2B
renege_abort() %>%
timeout(dist_screening)
branch_booking_none2B <- trajectory("no appointment needed") %>%
release("receptionist") %>%
set_attribute("booked_appt", 0) %>%
set_attribute("call_end", function() now(env_part2B)) # changed to env_part2B
branch_booking_phone2B <- trajectory("booking a f2f") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 1) %>%
set_attribute("call_end", function() now(env_part2B)) # changed to env_part2B
branch_booking_f2f2B <- trajectory("booking a f2f") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 2) %>%
set_attribute("call_end", function() now(env_part2B)) # changed to env_part2B
patient_flow2B <- trajectory("patient flow2B") %>%
branch(dist_contact_mode, TRUE,
branch_ecr,
branch_phone2B) %>%
timeout(dist_decision) %>%
branch(dist_booking_type, TRUE,
branch_booking_none2B,
branch_booking_phone2B,
branch_booking_f2f2B)
Populate environment using the trajectory for patient flow 2B
env_part2B %>%
add_resource("receptionist", capacity = sch_receptionists) %>%
add_generator("patient[8am-10am]", patient_flow2B,
from_to(t_hour(8), t_hour(10), dist_patient_interarrival1, every = t_day(1)),
mon = 2) %>%
add_generator("patient[10am-4pm]", patient_flow2B,
from_to(t_hour(10), t_hour(16), dist_patient_interarrival2, every = t_day(1)),
mon = 2) %>%
add_generator("patient[4pm-6pm]", patient_flow2B,
from_to(t_hour(16), t_hour(18), dist_patient_interarrival3, every = t_day(1)),
mon = 2)
## simmer environment: part 2B | now: 0 | next: 0
## { Monitor: in memory }
## { Resource: receptionist | monitored: TRUE | server status: 0(0) | queue status: 0(Inf) }
## { Source: patient[8am-10am] | monitored: 2 | n_generated: 0 }
## { Source: patient[10am-4pm] | monitored: 2 | n_generated: 0 }
## { Source: patient[4pm-6pm] | monitored: 2 | n_generated: 0 }
Run your model for 5 days
What is the average wait for calls to be answered, number abandoning their calls (hanging up) and the number booked for a face to face appointment in this simulated week?
env_part2B %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
summarise(countreneged = n())
| countreneged |
|---|
| 582 |
n=582
How many people get a F2F booking?
env_part2B %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
summarise(countF2F = n())
| countF2F |
|---|
| 160 |
n=160
What are the call waits?
part2B_callwait_results <- env_part2B %>%
get_mon_attributes() %>%
pivot_wider(-c(time, replication), names_from = key, values_from = value) %>%
filter(via_ecr == 0, arrival >= 0) %>%
mutate(wait_time_call = case_when(
!is.na(reneged) ~ 5,
is.na(call_taken) ~ (10 * 60) - arrival, # in queue/progress still, model run for 10 hours
TRUE ~ call_taken - arrival # has been answered
))
part2B_callwait_summary <- part2B_callwait_results %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(is.na(call_end) & is.na(reneged)),
countfinished = sum(!is.na(call_end) | !is.na(reneged)))
part2B_callwait_summary
| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|
| 3.825883 | 3.183036 | 3.479728 | 0 | 5 | 1.464796 | 5 | 0 | 1726 |
# # A tibble: 1 x 9
# median_wait mean_wait var_wait min_wait max_wait Qtile25th_wait Qtile75th_wait
#1 3.83 3.18 3.48 0 5 1.46 5
# countunfinished countfinished
# 0 1726
Now we are going to run our model 50 times and analyse different scenarios.
Use lapply to loop through and replicate the environment you created in Part 2 100 times
Update some of the code you used to look at results in pasrt 1 and 2. This time group_by(replication) before you summarise. E.g count of the reneged, then ungroup(). This will give the result for each replication. It is these that are then summarised again to give your final results for i.e. the mean of the mean number reneged in each replication.
envs_baseline %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
group_by(replication) %>%
summarise(countreneged = n()) %>%
ungroup() %>%
summarise(countreneged_mean = mean(countreneged))
| countreneged_mean |
|---|
| 579.86 |
n=580
envs_baseline %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
group_by(replication) %>%
summarise(countF2F = n()) %>%
ungroup() %>%
summarise(countF2F_mean = mean(countF2F))
| countF2F_mean |
|---|
| 143.04 |
n=143
Before we look at call waits let’s understand a little more about unfinished. using the staff schedule means that the final staff member finishes at 6pm on the 5th day so at (42460)+(18*60)= 6840, there are 4 in the system, 1 patient in progress and 3 waiting the staff member finishes that one at 6845 and then leaves with 3 remaining unfinished a similar thing occurs at the end of day 3:
get_mon_resources(envs_baseline) %>%
filter(replication == 1) %>%
tail(n = 10) # show last 10 rows
| resource | time | server | queue | capacity | queue_size | system | limit | replication | |
|---|---|---|---|---|---|---|---|---|---|
| 4184 | receptionist | 6824.155 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 4185 | receptionist | 6827.207 | 1 | 5 | 1 | Inf | 6 | Inf | 1 |
| 4186 | receptionist | 6829.155 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 4187 | receptionist | 6829.959 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| 4188 | receptionist | 6832.803 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 4189 | receptionist | 6835.244 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| 4190 | receptionist | 6835.407 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 4191 | receptionist | 6839.038 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| 4192 | receptionist | 6840.000 | 1 | 3 | 0 | Inf | 4 | Inf | 1 |
| 4193 | receptionist | 6845.235 | 0 | 3 | 0 | Inf | 3 | Inf | 1 |
get_mon_resources(envs_baseline) %>%
filter(replication == 1, time < 6240) %>%
tail(n = 10)
| resource | time | server | queue | capacity | queue_size | system | limit | replication | |
|---|---|---|---|---|---|---|---|---|---|
| 3326 | receptionist | 5387.154 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 3327 | receptionist | 5387.827 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| 3328 | receptionist | 5391.761 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 3329 | receptionist | 5392.154 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| 3330 | receptionist | 5392.374 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 3331 | receptionist | 5392.431 | 1 | 5 | 1 | Inf | 6 | Inf | 1 |
| 3332 | receptionist | 5393.999 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 3333 | receptionist | 5397.431 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| 3334 | receptionist | 5400.000 | 1 | 3 | 0 | Inf | 4 | Inf | 1 |
| 3335 | receptionist | 5401.615 | 0 | 3 | 0 | Inf | 3 | Inf | 1 |
envs_baseline %>%
get_mon_arrivals(ongoing = TRUE) %>% # 101972 rows
filter(replication == 1, is.na(end_time))
| name | start_time | end_time | activity_time | finished | replication |
|---|---|---|---|---|---|
| patient[4pm-6pm]212 | -1.000 | NA | NA | FALSE | 1 |
| patient[4pm-6pm]201 | 6790.286 | NA | NA | FALSE | 1 |
| patient[4pm-6pm]191 | 6774.405 | NA | NA | FALSE | 1 |
| patient[4pm-6pm]192 | 6775.585 | NA | NA | FALSE | 1 |
| patient[10am-4pm]1218 | -1.000 | NA | NA | FALSE | 1 |
| patient[8am-10am]653 | -1.000 | NA | NA | FALSE | 1 |
baseline_attributes <- envs_baseline %>%
get_mon_attributes() %>% #202732 obs
dplyr::select(-time) %>%
pivot_wider(names_from = "key", values_from = "value", values_fill = 0)
#101468 obs
baseline_callwait_results <- envs_baseline %>%
get_mon_arrivals(ongoing = TRUE) %>%
distinct() %>% # goes down to 101618, removes duplicates of unfinished
# still more in arrivals than attributes as attributes only has finished
full_join(baseline_attributes, by = c("name", "replication")) %>%
filter(via_ecr == 0) %>% # phone only also removes those yet to start e.g. start_time -1
mutate(wait_time_call = case_when(
# the caller reneged
reneged == 1 ~ 5,
# unfinished still in the system,
!finished ~ t_day(5) - start_time,
# everything else, (assume has been answered)
TRUE ~ pmax(end_time - start_time - activity_time, 0)
))
## calculate summary statistics
baseline_callwait_replication_summary <- baseline_callwait_results %>%
group_by(replication) %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(!finished),
countfinished = sum(finished)) %>%
mutate(scenario = "baseline")
# show top 10 rows
head(baseline_callwait_replication_summary, n = 10)
| replication | median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished | scenario |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 3.825883 | 3.183036 | 3.479728 | 0 | 5 | 1.464796 | 5 | 0 | 1726 | baseline |
| 2 | 3.744117 | 3.221769 | 3.258414 | 0 | 5 | 1.555424 | 5 | 0 | 1770 | baseline |
| 3 | 3.775021 | 3.116034 | 3.682359 | 0 | 5 | 1.206570 | 5 | 0 | 1772 | baseline |
| 4 | 3.952666 | 3.235613 | 3.446565 | 0 | 5 | 1.535425 | 5 | 0 | 1722 | baseline |
| 5 | 3.810451 | 3.201863 | 3.433046 | 0 | 5 | 1.446567 | 5 | 0 | 1736 | baseline |
| 6 | 3.550199 | 3.053321 | 3.688564 | 0 | 5 | 1.141632 | 5 | 0 | 1699 | baseline |
| 7 | 3.525437 | 3.005393 | 3.874786 | 0 | 5 | 1.005038 | 5 | 0 | 1652 | baseline |
| 8 | 3.808467 | 3.141373 | 3.616938 | 0 | 5 | 1.237069 | 5 | 0 | 1733 | baseline |
| 9 | 3.892609 | 3.168985 | 3.635974 | 0 | 5 | 1.225592 | 5 | 0 | 1750 | baseline |
| 10 | 3.623343 | 3.068928 | 3.714392 | 0 | 5 | 1.107349 | 5 | 0 | 1662 | baseline |
baseline_callwait_summary <- baseline_callwait_replication_summary %>%
summarise(median_wait = mean(median_wait),
mean_wait = mean(mean_wait),
var_wait = mean(var_wait),
min_wait = mean(min_wait),
max_wait = mean(max_wait),
Qtile25th_wait = mean(Qtile25th_wait),
Qtile75th_wait = mean(Qtile75th_wait),
countunfinished = mean(countunfinished),
countfinished = mean(countfinished)) %>%
mutate(scenario = "baseline") # to identify it later
baseline_callwait_summary
| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished | scenario |
|---|---|---|---|---|---|---|---|---|---|
| 3.765179 | 3.147497 | 3.580269 | 0 | 5 | 1.334405 | 5 | 0 | 1722.8 | baseline |
If a neighbouring practice closed leading to an increase in demand of 45%. Whilst the practice were able to better advertise the ECR route to patients such that 40% of all the patients (old and new) were able to use this route. What impact would this have on the call wait and reneging of patients?
dist_contact_mode_S <- function() sample(1:2, 1, FALSE, c(0.4, 0.6))
patient_flow3 <- trajectory("patient flow3") %>%
branch(dist_contact_mode_S, TRUE, # <- changed to the scenario distribution
branch_ecr,
branch_phone2) %>%
timeout(dist_decision) %>%
branch(dist_booking_type, TRUE,
branch_booking_none,
branch_booking_phone,
branch_booking_f2f)
dist_patient_interarrival1_S <- function() rexp(1, (120 * 1.45) / t_hour(2))
dist_patient_interarrival2_S <- function() rexp(1, (240 * 1.45) / t_hour(6))
dist_patient_interarrival3_S <- function() rexp(1, (40 * 1.45) / t_hour(2))
Repurpose some of the code you used to look at results in parts 1 and 2 of the workshop. Remember to store the results with a different name to your baseline results.
envs_scenario %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
group_by(replication) %>%
summarise(countreneged = n()) %>%
ungroup() %>%
summarise(countreneged_mean = mean(countreneged))
| countreneged_mean |
|---|
| 603.16 |
n=752
envs_scenario %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
group_by(replication) %>%
summarise(countF2F = n()) %>%
ungroup() %>%
summarise(countF2F_mean = mean(countF2F))
| countF2F_mean |
|---|
| 171.56 |
n=143
scenario_attributes <- envs_scenario %>%
get_mon_attributes() %>%
dplyr::select(-time) %>%
pivot_wider(names_from = "key", values_from = "value", values_fill = 0)
scenario_callwait_results <- envs_scenario %>%
get_mon_arrivals(ongoing = TRUE) %>%
distinct() %>%
full_join(scenario_attributes, by = c("name", "replication")) %>%
filter(via_ecr == 0) %>%
mutate(wait_time_call = case_when(
# the caller reneged
reneged == 1 ~ 5,
# unfinished still in the system,
!finished ~ t_day(5) - start_time,
# everything else, (assume has been answered)
TRUE ~ pmax(end_time - start_time - activity_time, 0)
))
## calculate summary statistics
scenario_callwait_replication_summary <- scenario_callwait_results %>%
group_by(replication) %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(!finished),
countfinished = sum(finished)) %>%
mutate(scenario = "scenario")
scenario_callwait_summary <- scenario_callwait_replication_summary %>%
summarise(median_wait = mean(median_wait),
mean_wait = mean(mean_wait),
var_wait = mean(var_wait),
min_wait = mean(min_wait),
max_wait = mean(max_wait),
Qtile25th_wait = mean(Qtile25th_wait),
Qtile75th_wait = mean(Qtile75th_wait),
countunfinished = mean(countunfinished),
countfinished = mean(countfinished)) %>%
mutate(scenario = "scenario")
# show top 10 results
head(scenario_callwait_summary, n = 10)
| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished | scenario |
|---|---|---|---|---|---|---|---|---|---|
| 3.890618 | 3.289315 | 3.122619 | 0.001257 | 5 | 1.598594 | 5 | 0 | 1747.56 | scenario |
Use t.test() function to compare waiting times between baseline and scenario and see if there is a significant difference in the results.
t.test(scenario_callwait_results$wait_time_call, baseline_callwait_results$wait_time_call)
##
## Welch Two Sample t-test
##
## data: scenario_callwait_results$wait_time_call and baseline_callwait_results$wait_time_call
## t = 16.071, df = 172339, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.1241180 0.1585979
## sample estimates:
## mean of x mean of y
## 3.290622 3.149264
This could actually be a paired t-test because they use the same common random number set and have the same number of replications therefore replication 1 of baseline and scenario are correlated (unless the scenario itself changes the andom number sampling). we should check the variances of the differences is less than the individual variances before using paired t test. To do it, we need to:
part3_summary_byreplication <- bind_rows(
baseline_callwait_replication_summary,
scenario_callwait_replication_summary
)
d <- with(
part3_summary_byreplication,
mean_wait[scenario == "baseline"] - mean_wait[scenario == "scenario"]
)
shapiro.test(d) # => p-value = 0.6375
##
## Shapiro-Wilk normality test
##
## data: d
## W = 0.95276, p-value = 0.04435
res <- t.test(mean_wait ~ scenario, data = part3_summary_byreplication, paired = TRUE)
res
##
## Paired t-test
##
## data: mean_wait by scenario
## t = -8.5021, df = 49, p-value = 3.295e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1753380 -0.1082975
## sample estimates:
## mean of the differences
## -0.1418177
p-value is greater than the significance level 0.05 implying that the distribution of the differences (d) are not significantly different from normal distribution. therefore we do a paired T-Test
data: mean_wait by scenario
t = -26.092, df = 49, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.4228384 -0.3623634
sample estimates:
mean of the differences
-0.3926009
The p-value of the test is < 2.2e-16, which is less than the significance level alpha = 0.05. We can then reject null hypothesis and conclude that the average wait in the baseline is significantly different from the average wait in the scenario
Using the ggplot2 package and geom_boxplot(), create boxplots of the means of the replications of the call waiting time in the baseline compared to the scenario. Set theme and title if you wish. To compare boxplots, you can either
a <- ggplot(baseline_callwait_replication_summary) +
geom_boxplot(aes(x = scenario, y = mean_wait),
outlier.shape = 19,
width = 0.3,
fill = "white",
colour = "orange") +
labs(y = "mean call wait (mins)") +
ylim(0, 4) +
theme(axis.title.x = element_blank(),
axis.text = element_text(size = 11))
b <- ggplot(scenario_callwait_replication_summary) +
geom_boxplot(aes(x = scenario, y = mean_wait),
outlier.shape = 19,
width = 0.3,
fill = "white",
colour = "orange") +
labs(y = "mean call wait (mins)") +
ylim(0, 4) +
theme(axis.title.x = element_blank(),
axis.text = element_text(size = 11))
grid.arrange(a, b, ncol = 2, nrow = 1, top = "Call wait mean results of baseline(left) and scenario(right)")
part3_summary <- bind_rows(baseline_callwait_summary, scenario_callwait_summary)
# plot
g <- ggplot() +
geom_boxplot(data = part3_summary_byreplication,
aes(x = scenario, y = mean_wait),
outlier.shape = 19,
width = 0.3,
fill = "white",
colour = "orange") +
geom_point(data = part3_summary,
aes(x = scenario, y = mean_wait),
shape = 19,
size = 2,
colour = "red") +
labs(x = "scenario", y = "mean call wait (mins)", title = "Boxplot of mean call wait from 50 replications") +
theme(axis.text = element_text(size = 11))
g
Looking at the results of baseline of 50 replications, what can you say about the results?
What do you interpret from the results of the scenario in comparison to the baseline?
What other analysis might you want to do of the model results?
Additional analysis: example if you want to compare ECR to Phones waits might do it this way:
baseline_wait_results <- envs_baseline %>%
get_mon_arrivals(ongoing = TRUE) %>%
distinct() %>%
full_join(baseline_attributes, by = c("name", "replication")) %>%
filter(start_time >= 0) %>% # so use this to remove the not started e.g. start_time -1
mutate(wait_time_call = case_when(
reneged == 1 ~ 5,
!finished ~ t_day(5) - start_time,
TRUE ~ pmax(end_time - start_time - activity_time, 0)
))
baseline_wait_replication_summary <- baseline_wait_results %>%
group_by(replication, via_ecr) %>% #group by both of these to get results for phone and ECR
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(!finished),
countfinished = sum(finished)) %>%
mutate(scenario = "baseline")
## `summarise()` has grouped output by 'replication'. You can override using the `.groups` argument.
# show top 10 results
head(baseline_wait_replication_summary, n = 10)
| replication | via_ecr | median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished | scenario |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 0 | 3.825883 | 3.183036 | 3.479728 | 0 | 5.000 | 1.464796 | 5.00000 | 0 | 1726 | baseline |
| 1 | 1 | 56.571437 | 121.432726 | 60680.123013 | 0 | 1129.364 | 14.716072 | 99.92010 | 3 | 354 | baseline |
| 2 | 0 | 3.744117 | 3.221769 | 3.258414 | 0 | 5.000 | 1.555424 | 5.00000 | 0 | 1770 | baseline |
| 2 | 1 | 65.218754 | 149.913667 | 69329.798636 | 0 | 1080.604 | 15.965540 | 149.12634 | 1 | 291 | baseline |
| 3 | 0 | 3.775021 | 3.116034 | 3.682359 | 0 | 5.000 | 1.206570 | 5.00000 | 0 | 1772 | baseline |
| 3 | 1 | 32.526169 | 64.839190 | 22473.421606 | 0 | 1075.412 | 5.002446 | 64.15030 | 6 | 314 | baseline |
| 4 | 0 | 3.952666 | 3.235613 | 3.446565 | 0 | 5.000 | 1.535425 | 5.00000 | 0 | 1722 | baseline |
| 4 | 1 | 41.584256 | 83.401127 | 35084.985697 | 0 | 1099.012 | 12.369562 | 73.55387 | 1 | 303 | baseline |
| 5 | 0 | 3.810451 | 3.201863 | 3.433046 | 0 | 5.000 | 1.446567 | 5.00000 | 0 | 1736 | baseline |
| 5 | 1 | 40.110651 | 76.963237 | 27211.750428 | 0 | 1008.028 | 9.223911 | 84.78485 | 0 | 299 | baseline |
baseline_wait_summary <- baseline_wait_replication_summary %>%
group_by(via_ecr) %>%
summarise(median_wait = mean(median_wait),
mean_wait = mean(mean_wait),
var_wait = mean(var_wait),
min_wait = mean(min_wait),
max_wait = mean(max_wait),
Qtile25th_wait = mean(Qtile25th_wait),
Qtile75th_wait = mean(Qtile75th_wait),
countunfinished = mean(countunfinished),
countfinished = mean(countfinished))
baseline_wait_summary
| via_ecr | median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 3.765179 | 3.147497 | 3.580269 | 0 | 5.000 | 1.334405 | 5.00000 | 0.00 | 1722.80 |
| 1 | 42.015194 | 94.007905 | 41412.497885 | 0 | 1052.733 | 10.070932 | 85.24734 | 4.08 | 302.48 |
What other ways might you want to present data?
For example, plot results from all replications
plot(get_mon_arrivals(envs_baseline), metric = "waiting_time")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 1062 rows containing non-finite values (stat_smooth).
## Warning: Removed 1062 row(s) containing missing values (geom_path).
plot(get_mon_arrivals(envs_scenario), metric = "waiting_time")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 29002 rows containing non-finite values (stat_smooth).
## Warning: Removed 29002 row(s) containing missing values (geom_path).
Histogram
ggplot(part3_summary_byreplication, aes(mean_wait)) +
geom_histogram(fill = "orange") +
labs(title = "Distribution of average call wait in 50 replications") +
facet_wrap(vars(scenario), ncol = 2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
What other opportunities are there for improvement of the process?
What other scenarios might you want to consider with the practice? What would your performance indicators be? How would you incorporate them in the model? Have a go if you have time. warm up period, staff rota, percentage coming via ecr, include preemption and restart
What else would need to be considered before they could implement changes suggested by the model? costs, other tasks the staff do, feasibility of suggested rota in real life
Are we running the model for a sufficient number of runs? process converted to R from: Simulations: The Practice of Model Development and Use, Robinson, S., 2nd Ed. 2014, Palgrave Macmillan
Remember back to our t test results we saved in res gave a confidence interval based on student’s t distribution
res$conf.int[1]
## [1] -0.175338
res$conf.int[2]
## [1] -0.1082975
Create a blank data frame to hold the cumulative results
check_baseline_replications <- data.frame(matrix(ncol = 5, nrow = 0))
x <- c("replication", "cmean", "tL", "tU", "percentdeviation")
colnames(check_baseline_replications) <- x
#for loop to see how the percent deviation of the confidence limits from the mean vary by the number of replications
for (i in 2:50) {
df <- baseline_callwait_replication_summary %>%
dplyr::select(replication, mean_wait) %>%
filter(replication <= i) %>%
summarise(replication = i,
cmean = mean(mean_wait),
tL = t.test(mean_wait)$conf.int[1],
tU = t.test(mean_wait)$conf.int[2],
percentdeviation = (cmean - tL) / cmean * 100) %>%
data.frame()
check_baseline_replications <- bind_rows(check_baseline_replications, df)
}
# show top 10 results
head(check_baseline_replications, n = 10)
| replication | cmean | tL | tU | percentdeviation |
|---|---|---|---|---|
| 2 | 3.202403 | 2.956327 | 3.448479 | 7.684104 |
| 3 | 3.173613 | 3.040728 | 3.306499 | 4.187198 |
| 4 | 3.189113 | 3.103887 | 3.274340 | 2.672417 |
| 5 | 3.191663 | 3.133636 | 3.249691 | 1.818099 |
| 6 | 3.168606 | 3.094869 | 3.242343 | 2.327118 |
| 7 | 3.145290 | 3.062986 | 3.227594 | 2.616748 |
| 8 | 3.144801 | 3.075910 | 3.213691 | 2.190611 |
| 9 | 3.147488 | 3.087915 | 3.207060 | 1.892702 |
| 10 | 3.139632 | 3.084423 | 3.194840 | 1.758443 |
| 11 | 3.153205 | 3.095464 | 3.210947 | 1.831200 |
To include the set_attribute using now() we reference the correct environment for each replication
We can now analyse results
envs_baselineB %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
group_by(replication) %>%
summarise(countreneged = n()) %>%
ungroup() %>%
summarise(countreneged_mean = mean(countreneged))
| countreneged_mean |
|---|
| 579.86 |
n=580
envs_baselineB %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
group_by(replication) %>%
summarise(countF2F = n()) %>%
ungroup() %>%
summarise(countF2F_mean = mean(countF2F))
| countF2F_mean |
|---|
| 143.04 |
n=143
baselineB_callwait_results <- envs_baselineB %>%
get_mon_attributes() %>%
pivot_wider(-c(time), names_from = key, values_from = value) %>%
filter(via_ecr == 0, arrival >= 0) %>%
mutate(wait_time_call = case_when(
!is.na(reneged) ~ 5,
is.na(call_taken) ~ (10 * 60) - arrival, # in queue/progress still, model run for 10 hours
TRUE ~ call_taken - arrival # has been answered
))
baselineB_callwait_replication_summary <- baselineB_callwait_results %>%
group_by(replication) %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(is.na(call_end) & is.na(reneged)),
countfinished = sum(!is.na(call_end) | !is.na(reneged)))
# show top 10 results
head(baselineB_callwait_replication_summary, n = 10)
| replication | median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 3.825883 | 3.183036 | 3.479728 | 0 | 5 | 1.464796 | 5 | 0 | 1726 |
| 2 | 3.744117 | 3.221769 | 3.258414 | 0 | 5 | 1.555424 | 5 | 0 | 1770 |
| 3 | 3.775021 | 3.116034 | 3.682359 | 0 | 5 | 1.206570 | 5 | 0 | 1772 |
| 4 | 3.952666 | 3.235613 | 3.446565 | 0 | 5 | 1.535425 | 5 | 0 | 1722 |
| 5 | 3.810451 | 3.201863 | 3.433046 | 0 | 5 | 1.446567 | 5 | 0 | 1736 |
| 6 | 3.550199 | 3.053321 | 3.688564 | 0 | 5 | 1.141632 | 5 | 0 | 1699 |
| 7 | 3.525437 | 3.005393 | 3.874786 | 0 | 5 | 1.005038 | 5 | 0 | 1652 |
| 8 | 3.808467 | 3.141373 | 3.616938 | 0 | 5 | 1.237069 | 5 | 0 | 1733 |
| 9 | 3.892609 | 3.168985 | 3.635974 | 0 | 5 | 1.225592 | 5 | 0 | 1750 |
| 10 | 3.623343 | 3.068928 | 3.714392 | 0 | 5 | 1.107349 | 5 | 0 | 1662 |
baselineB_callwait_summary <- baselineB_callwait_replication_summary %>%
summarise(median_wait = mean(median_wait),
mean_wait = mean(mean_wait),
var_wait = mean(var_wait),
min_wait = mean(min_wait),
max_wait = mean(max_wait),
Qtile25th_wait = mean(Qtile25th_wait),
Qtile75th_wait = mean(Qtile75th_wait),
countunfinished = mean(countunfinished),
countfinished = mean(countfinished)) %>%
mutate(scenario = "baseline") # to identify it later
baselineB_callwait_summary
| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished | scenario |
|---|---|---|---|---|---|---|---|---|---|
| 3.765179 | 3.147497 | 3.580269 | 0 | 5 | 1.334405 | 5 | 0 | 1722.8 | baseline |
baselineB_callwait_summary$mean_wait
## [1] 3.147497
baseline_callwait_summary$mean_wait
## [1] 3.147497
Trajectory must be included inside the lapply for replication
We can now analyse results
envs_scenarioB %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
group_by(replication) %>%
summarise(countreneged = n()) %>%
ungroup() %>%
summarise(countreneged_mean = mean(countreneged))
| countreneged_mean |
|---|
| 603.16 |
n=752
envs_scenarioB %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
group_by(replication) %>%
summarise(countF2F = n()) %>%
ungroup() %>%
summarise(countF2F_mean = mean(countF2F))
| countF2F_mean |
|---|
| 171.56 |
n=165
scenarioB_callwait_results <- envs_scenarioB %>%
get_mon_attributes() %>%
pivot_wider(-c(time), names_from = key, values_from = value) %>%
filter(via_ecr == 0, arrival >= 0) %>%
mutate(wait_time_call = case_when(
!is.na(reneged) ~ 5,
is.na(call_taken) ~ (10 * 60) - arrival, # in queue/progress still, model run for 10 hours
TRUE ~ call_taken - arrival # has been answered
))
scenarioB_callwait_replication_summary <- scenarioB_callwait_results %>%
group_by(replication) %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(is.na(call_end) & is.na(reneged)),
countfinished = sum(!is.na(call_end) | !is.na(reneged)))
# show top 10 results
head(scenarioB_callwait_replication_summary, n = 10)
| replication | median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 3.749845 | 3.209906 | 3.181770 | 0.0000000 | 5 | 1.464999 | 5 | 0 | 1666 |
| 2 | 3.932084 | 3.337168 | 3.053690 | 0.0000000 | 5 | 1.716950 | 5 | 0 | 1770 |
| 3 | 3.890046 | 3.282838 | 3.134920 | 0.0000000 | 5 | 1.540661 | 5 | 0 | 1789 |
| 4 | 3.764143 | 3.236900 | 3.134823 | 0.0000000 | 5 | 1.548852 | 5 | 0 | 1712 |
| 5 | 3.901525 | 3.279950 | 3.099232 | 0.0000000 | 5 | 1.651958 | 5 | 0 | 1723 |
| 6 | 4.060991 | 3.371249 | 3.040204 | 0.0111140 | 5 | 1.751085 | 5 | 0 | 1764 |
| 7 | 4.017161 | 3.352447 | 3.025200 | 0.0009612 | 5 | 1.774429 | 5 | 0 | 1763 |
| 8 | 3.904756 | 3.298098 | 3.105040 | 0.0000000 | 5 | 1.632506 | 5 | 0 | 1755 |
| 9 | 4.038012 | 3.353227 | 3.106330 | 0.0000000 | 5 | 1.716001 | 5 | 0 | 1764 |
| 10 | 3.684992 | 3.213686 | 3.210820 | 0.0000000 | 5 | 1.500257 | 5 | 0 | 1708 |
scenarioB_callwait_summary <- scenarioB_callwait_replication_summary %>%
summarise(median_wait = mean(median_wait),
mean_wait = mean(mean_wait),
var_wait = mean(var_wait),
min_wait = mean(min_wait),
max_wait = mean(max_wait),
Qtile25th_wait = mean(Qtile25th_wait),
Qtile75th_wait = mean(Qtile75th_wait),
countunfinished = mean(countunfinished),
countfinished = mean(countfinished)) %>%
mutate(scenario = "scenario") # to identify it later
scenarioB_callwait_summary
| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished | scenario |
|---|---|---|---|---|---|---|---|---|---|
| 3.890618 | 3.289315 | 3.122619 | 0.001257 | 5 | 1.598594 | 5 | 0 | 1747.56 | scenario |
scenarioB_callwait_summary$mean_wait
## [1] 3.289315
scenario_callwait_summary$mean_wait
## [1] 3.289315